﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace Bouyei.GeoCore.Geometries
{
    public class GeoBase
    {
        protected const double k0 = 1.57048761144159E-07;
        protected const double k1 = 5.05250178820567e-03;
        protected const double k2 = 2.98472900956587E-05;
        protected const double k3 = 2.14626669230084E-07;
        protected const double k4 = 2.22241238938534E-09;

        protected double e1 { get; private set; }
        protected double e2 { get; private set; }
        protected int precision { get; private set; }

        protected int zoneWide { get; private set; } = 3;//带宽
        protected double f { get; private set; }//扁率

        readonly double c = 0;
        const double degree = 57.295779513082320876798154814106;//180 / Math.PI;
        const double radian = 0.01745329251994329576923690768489;// Math.PI / 180;
        double m0 = 0;
        double m2 = 0;
        double m4 = 0;
        double m6 = 0;
        double m8 = 0;

        double a0 = 0;
        double a2 = 0;
        double a4 = 0;
        double a6 = 0;
        double a8 = 0;

        protected double lRadius { get; private set; }

        protected double sRadius { get; private set; }

        public GeoBase(double lRadius/*长半轴*/,
           double f/*递增变化量常量值(扁率)*/, int precision = 0)
        {
            this.precision = precision;
            this.lRadius = lRadius;
            this.f = f;
            this.sRadius = lRadius - lRadius * f;

            e1 = (Math.Pow(lRadius, 2) - Math.Pow(sRadius, 2)) / Math.Pow(lRadius, 2);
            e2 = (Math.Pow(lRadius, 2) - Math.Pow(sRadius, 2)) / Math.Pow(sRadius, 2);
            c = (lRadius * lRadius) / this.sRadius;

            m0 = lRadius * (1 - e1);
            m2 = 3 * e1 * m0 / 2.0;
            m4 = 5 * e1 * m2 / 4.0;
            m6 = 7 * e1 * m4 / 6.0;
            m8 = 9 * e1 * m6 / 8.0;

            a0 = m0 + m2 / 2.0 + 3 * m4 / 8.0 + 5 * m6 / 16.0 + 35 * m8 / 128.0;
            a2 = m2 / 2.0 + m4 / 2.0 + 15 * m6 / 32.0 + 7 * m8 / 16.0;
            a4 = m4 / 8.0 + 3 * m6 / 16.0 + 7 * m8 / 32.0;
            a6 = m6 / 32.0 + m8 / 16.0;
            a8 = m8 / 128.0;
        }

        public GeoBase(double lRadius = 6378137.0/*长半轴*/, double sRadius = 6356752.31414/*短半轴*/,
            int precision = 0, int zoneWide = 3)
        {
            this.precision = precision;
            this.lRadius = lRadius;
            this.sRadius = sRadius;
            this.f = 1 - sRadius / lRadius;
            this.zoneWide = zoneWide;

            e1 = (Math.Pow(lRadius, 2) - Math.Pow(sRadius, 2)) / Math.Pow(lRadius, 2);
            e2 = (Math.Pow(lRadius, 2) - Math.Pow(sRadius, 2)) / Math.Pow(sRadius, 2);
            c = (lRadius * lRadius) / this.sRadius;

            m0 = lRadius * (1 - e1);
            m2 = 3 * e1 * m0 / 2.0;
            m4 = 5 * e1 * m2 / 4.0;
            m6 = 7 * e1 * m4 / 6.0;
            m8 = 9 * e1 * m6 / 8.0;

            a0 = m0 + m2 / 2.0 + 3 * m4 / 8.0 + 5 * m6 / 16.0 + 35 * m8 / 128.0;
            a2 = m2 / 2.0 + m4 / 2.0 + 15 * m6 / 32.0 + 7 * m8 / 16.0;
            a4 = m4 / 8.0 + 3 * m6 / 16.0 + 7 * m8 / 32.0;
            a6 = m6 / 32.0 + m8 / 16.0;
            a8 = m8 / 128.0;

        }

        protected double ToDegree(double radian) => radian * degree;

        protected double ToRadian(double degree) => degree * radian;

        public virtual double Distance(Coordinate start, Coordinate end)
        {
            double L;
            bool isLB = start.X <= 360;// ((int)start.X / 1000000L) == 0;
            if (isLB)
            {
                //double sLat1 = ToRadian(start.Y);
                //double sLat2 = ToRadian(end.Y);
                //double a = sLat1 - sLat2;
                //double b = ToRadian(start.X - end.X);

                //s = 2 * Math.Asin(Math.Sqrt(Math.Pow(Math.Sin(a * 0.5), 2) +
                //        Math.Cos(sLat1) * Math.Cos(sLat2) * Math.Pow(Math.Sin(b * 0.5), 2)));
                //s *= lRadius;

                Coordinate s = new Coordinate()
                {
                    Y = ToRadian(start.Y),
                    X = ToRadian(start.X)
                };

                Coordinate d = new Coordinate()
                {
                    Y = ToRadian(end.Y),
                    X = ToRadian(end.X)
                };
                L = lRadius * Math.Acos(Math.Cos(s.Y) * Math.Cos(d.Y) * Math.Cos(s.X - d.X) + Math.Sin(s.Y) * Math.Sin(d.Y));

                return L;
            }
            else
            {
                L = Math.Sqrt(Math.Pow(start.X - end.X, 2) + Math.Pow(start.Y - end.Y, 2));
            }

            if (precision > 0)
                L = Math.Round(L, precision);

            return L;
        }

        /// <summary>
        /// 转换地理坐标系(高斯反算)
        /// </summary>
        /// <param name="geometry"></param>
        /// <returns></returns>
        public Geometry ToGeoCoordinates(Geometry geometry)
        {
            //检测第一个是否需要转换
            var coord = geometry.GetGemoetry(0)[0][0];
            if (coord.X > 0 && coord.X <= 360) return geometry;

            Geometry result = new Geometry(geometry.GeometryType);
            result.SRID = geometry.SRID;
            result.Geometric = new List<List<GeoSequence>>(geometry.GeometryCount);

            for (int i = 0; i < geometry.GeometryCount; ++i)
            {
                var geo = geometry.GetGemoetry(i);
                var r_geos = new List<GeoSequence>(geo.Count);

                for (int j = 0; j < geo.Count; ++j)
                {
                    var seq = geo[j];

                    var r_seq = new GeoSequence(seq.Count,
                        GeoOrientation.Clockwise, geometry.GeometryType);

                    for (int k = 0; k < seq.Count; ++k)
                    {
                        var rCoord = XYtoLB(seq[k]);
                        r_seq[k] = rCoord;
                    }
                    r_geos.Add(r_seq);
                }
                result.Geometric.Add(r_geos);
            }
            return result;
        }

        /// <summary>
        /// 转换平面坐标系(高斯正算)
        /// </summary>
        /// <param name="geometry"></param>
        /// <returns></returns>
        public Geometry ToPlaneCoordinates(Geometry geometry)
        {
            //检测第一个是否需要转换
            var coord = geometry.GetGemoetry(0)[0][0];
            if (coord.X < 0 || coord.X >= 10000000) return geometry;

            Geometry result = new Geometry(geometry.GeometryType);
            result.SRID = geometry.SRID;
            result.Geometric = new List<List<GeoSequence>>(geometry.GeometryCount);

            for (int i = 0; i < geometry.GeometryCount; ++i)
            {
                var o_geo = geometry.GetGemoetry(i);
                var r_geos = new List<GeoSequence>(o_geo.Count);

                for (int j = 0; j < o_geo.Count; ++j)
                {
                    var seq = o_geo[j];

                    var r_seq = new GeoSequence(seq.Count,
                        GeoOrientation.Clockwise, geometry.GeometryType);

                    for (int k = 0; k < seq.Count; ++k)
                    {
                        var rCoord = LBtoXY(seq[k]);
                        r_seq[k] = rCoord;
                    }
                    r_geos.Add(r_seq);
                }

                result.Geometric.Add(r_geos);
            }
            return result;
        }


        /// <summary>
        /// 高斯反算为大地经纬度,x经度，y纬度
        /// </summary>
        /// <param name="coord"></param>
        /// <returns></returns>
        public virtual Coordinate XYtoLB(Coordinate coord)
        {
            int proNo = 0;
            if (zoneWide == 3)
                proNo = (int)((coord.X - 1.5) / 1000000L);
            else
                proNo = (int)(coord.X / 1000000L);

            double centerl = 0;
            if (zoneWide == 3)
                centerl = proNo * zoneWide;
            else
                centerl = (proNo - 1) * zoneWide + (zoneWide / 2);

            centerl = ToRadian(centerl);

            double y1 = coord.X - 500000L - proNo * 1000000L;
            double e = k0 * coord.Y;

            double se = Math.Sin(e);
            double bf = e + Math.Cos(e) * (k1 * se - k2 * Math.Pow(se, 3) + k3 * Math.Pow(se, 5) - k4 * Math.Pow(se, 7));

            double t = Math.Tan(bf);
            double t2 = Math.Pow(t, 2);

            double n2 = e2 * Math.Pow(Math.Cos(bf), 2);
            double v2 = 1 + n2;
            double v = Math.Sqrt(v2);
            double N = c / v;
            double YdN = y1 / N;
            double _bfc = 1 / Math.Cos(bf);
            double t4 = Math.Pow(t2, 2);
            double YdN2 = Math.Pow(YdN, 2);
            double YdN4 = Math.Pow(YdN2, 2);

            double By = bf - (v2 * t) * YdN2 * 0.5
                + (5 + 3 * t2 + n2 - 9 * n2 * t2) * (v2 * t) * YdN4 / 24.0
                - (61 + 90 * t2 + 45 * t4) * (v2 * t) * YdN4 * YdN2 / 720.0;

            double Lx = _bfc * YdN
                - (1 + 2 * t2 + n2) * _bfc * YdN2 * YdN / 6.0
                + (5 + 28 * t2 + 24 * t4 + 6 * n2 + 8 * n2 * t2) * _bfc * YdN4 * YdN / 120.0
                + centerl;

            return new Coordinate()
            {
                X = ToDegree(Lx),
                Y = ToDegree(By)
            };
        }

        /// <summary>
        /// 高斯正算为XY平面坐标
        /// </summary>
        /// <param name="coord"></param>
        /// <returns></returns>
        //public virtual Coordinate LBtoXY(Coordinate coord)
        //{
        //    int ProjNo = (int)((coord.X + 1.5) / zoneWide); //6度带=(Y+6)/6
        //    double centerl = ProjNo * zoneWide;//6度带=ProjNo*zoneWid-3
        //    centerl = ToRadian(centerl);

        //    double x = ToRadian(coord.X); //经度转换为弧度
        //    double y = ToRadian(coord.Y); //纬度转换为弧度
        //    double ex = 2 * f - f * f;
        //    double ee = ex * (1.0 - ex);

        //    var latsin1 = Math.Sin(y);
        //    var lattan1 = Math.Tan(y);
        //    var latcos1 = Math.Cos(y);

        //    double NN = lRadius / Math.Sqrt(1.0 - ex * latsin1 * latsin1);
        //    double T = lattan1 * lattan1;
        //    double C = ee * latcos1 * latcos1;
        //    double A = (x - centerl) * latcos1;

        //    double M = lRadius * ((1 - ex / 4.0 - 3 * ex * ex / 64.0 - 5 * ex * ex * ex / 256.0) * y
        //        - (3 * ex / 8.0 + 3 * ex * ex / 32.0 + 45 * ex * ex * ex / 1024) * Math.Sin(2 * y)
        //     + (15 * ex * ex / 256.0 + 45 * ex * ex * ex / 1024.0) * Math.Sin(4 * y) - (35 * ex * ex * ex / 3072.0) * Math.Sin(6 * y));

        //    double xval = NN * (A + (1 - T + C) * A * A * A / 6 + (5 - 18 * T + T * T + 72 * C - 58 * ee) * A * A * A * A * A / 120.0);

        //    double yval = M + NN * Math.Tan(y) * (A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24.0
        //     + (61 - 58 * T + T * T + 600 * C - 330 * ee) * A * A * A * A * A * A / 720.0);
        //    double X0 = 1000000L * (ProjNo) + 500000L;

        //    return new Coordinate()
        //    {
        //        X = xval + X0,
        //        Y = yval
        //    };
        //}

        /// <summary>
        /// 高斯正算为XY平面坐标,提高精度
        /// </summary>
        /// <param name="coord"></param>
        /// <returns></returns>
        public virtual Coordinate LBtoXY(Coordinate coord)
        {
            int ProjNo = 0;
            if (zoneWide == 3)
                ProjNo = (int)((coord.X + 1.5) / zoneWide);
            else
                ProjNo = (int)(coord.X / zoneWide);

            double L0 = 0;
            if (zoneWide == 3)
                L0 = ProjNo * zoneWide;
            else
                L0 = ProjNo * zoneWide + (zoneWide / 2);

            double B = ToRadian(coord.Y);
            double L = ToRadian(coord.X - L0);
            double t = Math.Tan(B);
            double t2 = Math.Pow(t, 2);
            double L2 = Math.Pow(L, 2);

            double cosb = Math.Cos(B);
            double sinb = Math.Sin(B);

            double sinb2 = Math.Pow(sinb, 2);
            double cosb2 = Math.Pow(cosb, 2);
            double t4 = Math.Pow(t2, 2);
            double n2 = e2 * cosb2;

            double N = lRadius / Math.Sqrt(1 - e1 * sinb2);
            double Y = a0 * B - sinb * cosb * ((a2 - a4 + a6) + (2 * a4 - 16 * a6 / 3.0) * sinb2 + 16* a6 * Math.Pow(sinb2, 2) / 3.0 );

            double y = Y + N * t * cosb2 * L2 * (
                0.5 + cosb2 * L2 * ((5 - t2 + 9 * n2 + 4 * Math.Pow(n2, 2)) / 24.0 +cosb2 * L2 * (61 - 58 * t2 + t4) / 720.0)
                );

            double x = N * cosb * L * (
                1 + cosb2 * L2 * ((1 - t2 + n2) / 6.0+ (5 - 18 * t2 + t4 + 14 * n2 - 58 * n2 * t2) * cosb2 * L2 / 120.0)
            );

            x = x + 500000 + ProjNo * 1000000;

            return new Coordinate()
            {
                X = x,
                Y = y
            };
        }

        /// <summary>
        /// 距离插值
        /// </summary>
        /// <param name="coords"></param>
        /// <param name="distMeters"></param>
        public virtual List<Coordinate> InterpolationPoints(List<Coordinate> coords,double distMeters)
        {
            List<Coordinate> newList = new List<Coordinate>(coords.Count << 1);

            int len = coords.Count - 1;
            for (int i = 0; i < len; ++i)
            {
                var start = coords[i];
                var end = coords[i + 1];

                newList.Add(start);
     
                var dist = start.Distance(end);
                var cnt = (int)(dist / distMeters);
                if (cnt == 0)
                {
                    newList.Add(end);
                    continue;
                }
                cnt += 1;

                double segX = Math.Abs(end.X - start.X) / cnt;
                double segY = Math.Abs(end.Y - start.Y) / cnt;

                for (int j = 1; j < cnt; ++j)
                {
                    var center = new Coordinate(start.X + segX * j, start.Y + segY * j);
                    newList.Add(center);
                }
                newList.Add(end);
            }
            return newList;
        }

        /// <summary>
        /// 简化集合
        /// </summary>
        /// <param name="coords"></param>
        /// <param name="tolerance"></param>
        /// <returns></returns>
        public virtual List<Coordinate> SimplifyPoints(List<Coordinate> coords, double tolerance)
        {
            bool[] usedList = new bool[coords.Count];

            SimplifyTo(usedList, coords, 0, coords.Count - 1, tolerance);

            List<Coordinate> list = new List<Coordinate>(coords.Count);
            for (int i = 0; i < usedList.Length; ++i)
            {
                if (usedList[i])
                    list.Add(coords[i]);
            }
            return list;
        }

        private void SimplifyTo(bool[] usedList,List<Coordinate> Nodes,int sIndex,int eIndex,double tolerance)
        {
            if (eIndex == sIndex + 1) return;

            Coordinate from = Nodes[sIndex], to = Nodes[eIndex];

            //var sqrt = Math.Sqrt(Math.Pow(from.Y - to.Y, 2) + Math.Pow(from.X - to.X, 2));
            var x_dist = (to.X - from.X);//方向很重要
            var y_dist = (from.Y - to.Y);//方向很重要
            var xy_dist = ((from.X * to.Y) - (to.X * from.Y));

            var xysqrt = Math.Sqrt(Math.Pow(x_dist, 2) + Math.Pow(y_dist, 2));

            double max = 0;
            int maxIndex = sIndex;

            //计算点到指定线段的最大距离
            for (int i = sIndex + 1; i < eIndex; ++i)
            {
                var dist = PointToSegment(Nodes[i],x_dist, y_dist, xy_dist, xysqrt);
                if (dist > max)
                {
                    max = dist;
                    maxIndex = i;
                }
            }

            //删除该范围点
            if (max <= tolerance)
            {
                for (int i = sIndex + 1; i < eIndex; ++i)
                {
                    usedList[i] = true;
                }
            }
            else
            {
                SimplifyTo(usedList,Nodes,sIndex, maxIndex,tolerance);
                SimplifyTo(usedList,Nodes,maxIndex, eIndex,tolerance);
            }
        }
     
        private double PointToSegment(Coordinate coord,double x_dist, double y_dist, double xy_dist,
            double xysqrt)
        {
            var dist = Math.Abs(x_dist * coord.Y + y_dist * coord.X + xy_dist) / xysqrt;

            return dist;
        }
    }
}
